#ifndef GEO3D_H
#define GEO3D_H

#include "utils.h"

namespace geo3d {
    struct Point3d {
        double x, y, z;
        Point3d() {}
        Point3d(double x, double y, double z) : x(x), y(y), z(z) {}

        // 重载运算符：向量加法
        Point3d operator+(const Point3d& p) const {
            return Point3d(x + p.x, y + p.y, z + p.z);
        }

        // 重载运算符：向量减法
        Point3d operator-(const Point3d& p) const {
            return Point3d(x - p.x, y - p.y, z - p.z);
        }

        Point3d& operator+=(const Point3d& p) {
            x += p.x;
            y += p.y;
            z += p.z;
            return *this;
        }

        Point3d& operator-=(const Point3d& p) {
            x -= p.x;
            y -= p.y;
            z -= p.z;
            return *this;
        }

        // 重载运算符：向量与标量乘法
        Point3d operator*(double d) const {
            return Point3d(x * d, y * d, z * d);
        }
        Point3d operator/(double d) const {
            return Point3d(x / d, y / d, z / d);
        }
    };

    typedef struct Line {
        Point3d a;
        Point3d b;
        Line() {}
        Line(Point3d a, Point3d b): a(a), b(b) {}
    } Line;

    typedef struct Fault {
        std::vector<std::vector<Point3d>> lines;
        std::string name;
        INDEX horizontal_size; // 断层柱线段数
        INDEX vertical_size; // z方向点数
    } Fault;

    // 线性插值函数
    // 当 t=0 时返回 A，当 t=1 时返回 B，中间值返回 A 与 B 的线性组合
    Point3d interpolate(const Point3d &A, const Point3d &B, double t);

    // 对3d线段进行n等分
    void splitLine(Point3d &a, Point3d &b, int n, std::vector<Point3d> &points); 

    // 查找最接近点p的点索引
    void findClosestPoint(const std::vector<Point3d> &surface_points, const Point3d &p, INDEX &index, double &dist);

    // 查找最接近给出柱线的柱线索引
    void findClosestLine(const std::vector<Point3d> &bottom_points, const std::vector<Point3d> &top_points, const Point3d &bottom_p, const Point3d &top_p, INDEX &index, double &dist);

    // 利用指数方法在四面体内采样一个点
    Point3d samplePointInTetrahedron(const Point3d &A, const Point3d &B, const Point3d &C, const Point3d &D, std::mt19937 &gen);

    double distance(Point3d a, Point3d b);
    double distance2(Point3d a, Point3d b);
    double magnitude(Point3d p);
    Point3d crossProduct(Point3d a, Point3d b);
    double dotProduct(Point3d a, Point3d b);
    // ab, ac 夹角余弦值
    double cosAngle(Point3d a, Point3d b, Point3d c);
    double triangleArea(Point3d a, Point3d b, Point3d c);
    // Compute the incenter of the tetrahedron
    Point3d computeIncenter(Point3d a, Point3d b, Point3d c, Point3d d);
    // Compute the mass center of the tetrahedron
    Point3d computeMasscenter(Point3d a, Point3d b, Point3d c, Point3d d);
    // 直线 line 是否与 box 有交点
    bool lineBoxIntersection(const Line& line, const Point3d& min_bnd, const Point3d& max_bnd);
    // 直线 line 是否与三角形 abc 有交点，有则返回交点 P, 返回 -1 : 平行, 1 : 交点在三角形内部, 0 : 交点在三角形外部
    int lineTriangleIntersection(const Line& line, const Point3d& a, const Point3d& b, const Point3d& c, Point3d& P);
}; // namespace geo3d

#endif // GEO3D_H
